home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C/C++ Users Group Library 1996 July
/
C-C++ Users Group Library July 1996.iso
/
vol_200
/
275_02
/
prob22.c
< prev
next >
Wrap
Text File
|
1980-01-01
|
28KB
|
971 lines
/* prob22.c */
/* program for LCA22 option 't' */
/* calculate probabilities related to evolution */
/* Harold V. McIntosh, 10 August 1987 */
/* 29 February 1988 - adapted from LCA21 to LCA22 */
/* references: */
/* */
/* W. John Wilbur, David J. Lipman and Shihab A. Shamma */
/* On the prediction of local patterns in cellular automata */
/* Physica 19D 397-410 (1986) */
/* */
/* Howard A. Gutowitz, Jonathan D. Victor and Bruce W. Knight */
/* Local structure theory for cellular automata */
/* Physica 28D 18-48 (1987) */
/* Copyright (C) 1987 */
/* Copyright (C) 1988 */
/* Harold V. McIntosh */
/* Gerardo Cisneros S. */
# define BROW 3 /* row for bar charts */
# define EROW 1 /* row for evolution synopsis */
# define AORG 0 /* x-origin of contour plot */
# define BORG 109 /* x-origin of 2-block param */
# define CORG 219 /* x-origin of Bernstein plot */
/* edit the probability screen */
edtri() {char c;
videomode(COLGRAF);
videopalette(YELREGR);
while (0<1) {
videocursor(0,0,0);
videoputc(':',2);
scrrul();
videocursor(0,0,36);
videoputc('?',2);
c=kbdin();
if (c == '\015') break;
videocursor(0,0,38);
videoputc(c,2);
videocursor(0,0,36);
videoputc(' ',0);
switch (c) {
case '+': videopalette(WHCYMAG); break;
case '-': videopalette(YELREGR); break;
case 'a': asfreq(3); break;
case 'e': pevolve(); break;
case 'g': lifreq(50,2); break;
case 'G': lifreq(200,1); break;
case 'm': moncar(1,2); break;
case 'M': moncar(50,1); break;
case 'x': pdiff(100); break;
case 'i': pidiff(100); break;
case 'j': pjdiff(100); break;
case 'y': pdiff3(100); break;
case 'z': pdiff4(100); break;
case 't': twoblockfreq(100); break;
case '1': nblclr(); oneblfreq(8*BROW,300,48); break;
case '2': nblclr(); twoblfreq(8*BROW,300,48); break;
case '3': nblclr(); thrblfreq(8*BROW,300,48); break;
case '4': nblclr(); foublfreq(8*BROW,300,48); break;
case '5': nblclr(); fivblfreq(8*BROW,300,48); break;
case '6': nblclr(); sixblfreq(8*BROW,300,48); break;
case '?': videomode(COLGRAF); videopalette(YELREGR); trmenu(); break;
case '/': videomode(COLGRAF); videopalette(YELREGR); break;
default: break;
}; /* end switch */
}; /* end while */
videopalette(WHCYMAG);
videomode(T80X25);
} /* end edtri */
/* show menu */
trmenu() {
videoscroll(BROW,0,BROW+8,40,0,0);
videocursor(0,BROW,0);
printf("a - a priori estimates\n");
printf("m,M,g,G - sample evolution\n");
printf("xyz - selfconsistent probabilities\n");
printf("xij - iterated s-c probabilities\n");
printf("t - graph 2 block probabilities\n");
printf("123456 - n-block bar charts\n");
printf("+- - change color pallette\n");
printf("e - 12 lines evolution\n");
printf("/?(clear screen, show menu), <cr>(exit)\n");
}
/* show twelve lines of evolution at top of screen */
pevolve() {int i, j;
videoscroll(EROW,0,EROW+1,40,0,0);
asctobin();
for (j=8*EROW; j<8*(EROW+2)-3; j++) {
for (i=0; i<AL; i++) videodot(j,i,arr1[i]);
onegen(AL);
};
}
/* Clear a space for the n-block frequencies */
nblclr() {videoscroll(BROW,0,BROW+8,40,0,0);}
/* make a frame for a graph */
/* (x,y) = lower left corner; e.g. (0,0) */
/* n = dimension of frame */
gfram(x,y,n) int x, y, n; {int i;
for (i=0; i<=n; i++) videodot(199-y-i,x,0);
for (i=0; i<=n; i++) videodot(199-y-i,x+n,0);
for (i=0; i<=n; i++) videodot(199-n-y,x+i,0);
for (i=0; i<=n; i++) videodot(199-y,x+i,0);
for (i=0; i<=n; i+=10) videodot(199-y-i,x,3);
for (i=0; i<=n; i+=10) videodot(199-y-i,x+n,3);
for (i=0; i<=n; i+=10) videodot(199-n-y,x+i,3);
for (i=0; i<=n; i+=10) videodot(199-y,x+i,3);
}
/* put a diagonal in a graph */
gdiag(x,y,n) int x, y, n; {int i;
for (i=0; i<=n; i+=2) videodot(199-y-i,x+i,3);
}
/* graph Bernstein polynomial */
bgraf(x,y,k,n) int x, y, k, n; {int i; double bern(), en, dp, p;
if (n==0) return;
en=(double)(n); dp=1.0/en;
for (i=0,p=0.0; i<n; i++,p+=dp) {videodot(199-y-(int)(en*bern(p,k)),x+i,1);};
}
/* "Monte Carlo" determination of probabilities */
moncar(n,l) int n, l; {
int i, j, k, b[KK], bb[KK][KK];
double bf[KK], bbf[KK][KK];
nblclr();
gfram(BORG,0,100);
asctobin();
for (k=0; k<n; k++) {
onegen(AL);
for (i=0; i<KK; i++) b[i]=0;
for (i=0; i<AL; i++) b[arr1[i]]++;
for (i=0; i<KK; i++) bf[i]=((double)(b[i]))/((double)(AL));
for (i=0; i<KK; i++) for (j=0; j<KK; j++) bb[i][j]=0;
for (i=1; i<AL; i++) bb[arr1[i-1]][arr1[i]]++;
bb[arr1[AL-1]][arr1[0]]++;
for (i=0; i<KK; i++) for (j=0; j<KK; j++)
bbf[i][j]=((double)(bb[i][j]))/((double)(AL));
videodot(199-(int)(100.0*bbf[1][1]),BORG+(int)(100.0*bf[1]),l);
};
videocursor(0,BROW+7,0);
printf("(Monte Carlo) ");
for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,bf[i]);
videocursor(0,BROW+8,0);
for (i=0; i<KK; i++) for (j=0; j<KK; j++)
printf("%1d%1d:%5.3f ",i,j,bbf[i][j]);
}
/* Generate coefficients of 2nd generation Bernstein Polynomial */
berncoef() {
int i, i0, i1, i2, i3, i4;
for (i=0; i<BD; i++) bp[i]=0.0;
for (i0=0; i0<KK; i0++) {
for (i1=0; i1<KK; i1++) {
for (i2=0; i2<KK; i2++) {
for (i3=0; i3<KK; i3++) {
for (i4=0; i4<KK; i4++) {
if (ascrule[i0][i1][i2][i3][i4]=='1') bp[i0+i1+i2+i3+i4]+=1.0;
};};};};};
}
/* Generate coefficients of 3rd generation Bernstein Polynomial */
bernthrd() {
int i, i0, i1, i2, i3, i4, i5, i6, i7, i8;
int j0, j1, j2, j3, j4;
asctobin();
for (i=0; i<BD; i++) bp[i]=0.0;
for (i0=0; i0<KK; i0++) {
for (i1=0; i1<KK; i1++) {
for (i2=0; i2<KK; i2++) {
for (i3=0; i3<KK; i3++) {
for (i4=0; i4<KK; i4++) {
for (i5=0; i5<KK; i5++) {
for (i6=0; i6<KK; i6++) {
for (i7=0; i7<KK; i7++) {
for (i8=0; i8<KK; i8++) {
j0=binrule[i0][i1][i2][i3][i4];
j1=binrule[i1][i2][i3][i4][i5];
j2=binrule[i2][i3][i4][i5][i6];
j3=binrule[i3][i4][i5][i6][i7];
j4=binrule[i4][i5][i6][i7][i8];
if (ascrule[j0][j1][j2][j3][j4]=='1') bp[i0+i1+i2+i3+i4+i5+i6+i7+i8]+=1.0;
};};};};};};};};};
}
/* Generate coefficients of 4th generation Bernstein Polynomial */
bernfrth() {
int i, i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, ia, ib, ic;
int j0, j1, j2, j3, j4, j5, j6, j7, j8;
int k0, k1, k2, k3, k4;
asctobin();
for (i=0; i<BD; i++) bp[i]=0.0;
for (i0=0; i0<KK; i0++) {
for (i1=0; i1<KK; i1++) {
for (i2=0; i2<KK; i2++) {
for (i3=0; i3<KK; i3++) {
for (i4=0; i4<KK; i4++) {
for (i5=0; i5<KK; i5++) {
for (i6=0; i6<KK; i6++) {
for (i7=0; i7<KK; i7++) {
for (i8=0; i8<KK; i8++) {
for (i9=0; i9<KK; i9++) {
for (ia=0; ia<KK; ia++) {
for (ib=0; ib<KK; ib++) {
for (ic=0; ic<KK; ic++) {
j0=binrule[i0][i1][i2][i3][i4];
j1=binrule[i1][i2][i3][i4][i5];
j2=binrule[i2][i3][i4][i5][i6];
j3=binrule[i3][i4][i5][i6][i7];
j4=binrule[i4][i5][i6][i7][i8];
j5=binrule[i5][i6][i7][i8][i9];
j6=binrule[i6][i7][i8][i9][ia];
j7=binrule[i7][i8][i9][ia][ib];
j8=binrule[i8][i9][ia][ib][ic];
k0=binrule[j0][j1][j2][j3][j4];
k1=binrule[j1][j2][j3][j4][j5];
k2=binrule[j2][j3][j4][j5][j6];
k3=binrule[j3][j4][j5][j6][j7];
k4=binrule[j4][j5][j6][j7][j8];
i=i0+i1+i2+i3+i4+i5+i6+i7+i8+i9+ia+ib+ic;
if (ascrule[k0][k1][k2][k3][k4]=='1') bp[i]+=1.0;
};};};};};};};};};};};};};
}
/* evaluate the nth generation Bernstein polynomial at point p */
double bern(p,n) double p; int n; {double q, s, x, r; int i, d;
d=4*n-3;
if (p>0.99) return bp[d];
q=1.0-p; r=p/q; s=0.0; x=1.0;
for (i=0; i<d; i++) x*=q;
for (i=0; i<=d; i++, x*=r) s+=bp[i]*x;
s*=0.9998;
return(s+0.0001);
}
/* graph the probability of the second generation */
pdiff(n) int n; {
berncoef();
gfram(CORG,0,n);
gdiag(CORG,0,n);
bgraf(CORG,0,2,n);
}
/* graph the iterated probability of the second generation */
pidiff(n) int n; {int i; double bern(), en, p;
en=1.0/((double)(n));
berncoef();
gdiag(CORG,0,n);
for (i=0; i<=n; i++) {
p=((double)(i))*en;
videodot(199-(int)(100.0*bern(bern(p,